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^ ■ We consider the formation and evolution of vortices in a hydrodynamic 

<^ ■ shearing-sheet model. The evolution is done numerically using a version of the 

ZEUS code. Consistent with earlier results, an injected vorticity field evolves 
into a set of long-lived vortices each of which has radial extent comparable to 
^ ■ the local scale height. But we also find that the resulting velocity field has pos- 

itive shear stress {T,5vr5v^). This effect appears only at high resolution. The 
O ' transport, which decays with time as t~^/^, arises primarily because the vortices 

Q ■ drive compressive motions. This result suggests a possible mechanism for angu- 

^ ■ lar momentum transport in low-ionization disks, with two important caveats: a 

mechanism must be found to inject vorticity into the disk, and the vortices must 
"p i' not decay rapidly due to three-dimensional instabilities. 
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C^ ' 1. Introduction 

Astrophysical disks are common because the specific angular momentum of the mat- 
ter inside them is well-conserved. They evolve because angular momentum conservation is 
weakly compromised, either because of diffusion of angular momentum within the disk or 
because of direct application of external torques. 

In astrophysical disks composed of a well-ionized plasma it is likely that some, perhaps 
most, of the evolution is driven by diffusion of angular momentum within the disk. This view 
is certainly consistent with observations of steadily accreting cataclysmic variable systems 
like UX Ursa Majoris (Baptista et al. 1998; Baptista 2004), whose radial surface-brightness 
profile is consistent with steady accretion-flow models in which the bulk of the accretion 
energy is dissipated within the disk. 
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Angular momentum diffusion in well-ionized disks is likely driven by magnetohydrody- 
namic (MHD) turbulence. Analytic analyses, numerical experiments, and laboratory evi- 
dence strongly suggest that well-coupled plasmas in differentially-rotating ffows are subject 
to the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998; Balbus 2003). But 
MHD turbulence is initiated by the MRI only so long as the plasma is sufficiently ionized 
to couple to the magnetic field (Kunz & Balbus 2004; Desch 2004). In disks around young 
stars, cataclysmic-variable and X-ray binary disks in quiescence, and possibly the outer 
parts of AGN disks, the plasma may be too neutral to support magnetic activity (Gammie 
& Menou 1998; Menou 2000; Stone, Gammie, Balbus, & Hawley 2000; Menou & Quataert 
2001; Fromang et al. 2002). This motivates interest in non-MHD angular momentum trans- 
port mechanisms. 

Within the last few years, a body of work has been developed suggesting that vortices can 
be generated as a result of global hydrodynamic instability (Hawley 1987; Blaes & Hawley 
1988; Hawley 1990; Lovelace, Li, Colgate, & Nelson 1999; Li, Finn, Lovelace, & Colgate 
2000) or local hydrodynamic instability (Klahr & Bodenheimer 2003), that vortices in disks 
may be long-lived (Godon & Livio 1999, 2000; Umurhan & Regev 2004; Barranco & Marcus 
2005), and that these vortices may be related to an outward flux of angular momentum (Li, 
Colgate, Wendroff, & Liska 2001; Barranco & Marcus 2005). If these claims can be verifled 
then the consequences for low-ionization disks would be profound. 

Here we investigate the evolution of a disk that is given a large initial vortical velocity 
perturbation. Our study is done in the context of a (two-dimensional) shearing-sheet model, 
which permits us to resolve the dynamics to a degree that is not currently possible in a global 
disk model. Our model is also fully compressible, unlike previous work using a local model 
(Umurhan & Regev 2004; Barranco & Marcus 2005). The former assume incompressible 
flow and the latter use the anelastic approximation (e.g.. Cough 1969), which fllters out the 
high-frequency acoustic waves. We will show that compressibility and acoustic waves play 
an essential part in the angular momentum transport. 

Our paper is organized as follows. In §2 we describe the model. In §3 we describe 
the evolution of a flducial, high-resolution model. In §4 we investigate the dependence of 
the results on model parameters. And in §5 we describe implications, with an emphasis on 
key open questions: are the vortices destroyed by three-dimensional instabilities?; and do 
mechanisms exist that can inject vorticity into the disk? 
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2. Model 

The shearing-sheet model is obtained via a rigorous expansion of the two-dimensional 
hydrodynamic equations of motion to lowest order in H/R, where H = Cs/fi is the disk 
scale height (c^ is the isothermal sound speed and Q is the local rotation frequency) and R 
is the local radius. See Narayan et al. (1987) for a description. Adopting a local Cartesian 
coordinate system where the x axis is oriented parallel to the radius vector and the y axis 
points forward in azimuth, the equations of motion become 

- + EV.. = 0. (1) 

dv VP _ ^9 . 

— + — ^ + 2nxv- 2qn^xx = 0, (2) 

where S and P are the two-dimensional density and pressure, v is the fluid velocity and 
d/dt is the Lagrangian derivative. The third and fourth terms in equation (2) represent the 
Coriolis and centrifugal forces in the local model expansion, where q = —(1/2) dlnQ'^/dlnr is 
the shear parameter. We will assume throughout that q = 3/2, corresponding to a Keplerian 
shear profile. We close the above equations with an isothermal equation of state 

P = c% (3) 

where Cg is constant in time and space. 

Equations (1) through (3) can be combined to show that the vertical component of 

potential vorticity 

^ ^ {Vxv + 2n)-z ^^^ 

Zj 

is a constant of the motion; i.e., the potential vorticity of fluid elements in two dimensions 
is conserved. 

An equilibrium solution to the equations of motion is 

S = So = const. (5) 

P = c^Sq = const. (6) 

V. = (7) 

Vy = —qQx (8) 

Thus the differential rotation of the disk makes an appearance in the form of a linear shear. 

We integrate the above equations using a version of the ZEUS code (Stone & Norman 
1992). ZEUS is a time-explicit, operator-split scheme on a staggered mesh. It uses artificial 
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viscosity to capture shocks. Our computational domain is a rectangle of size Lx x Ly con- 
taining A^^ X Ny grid cells. The numerical resolution is therefore Ax x Ay = L^j^x x Ly/Ny. 

Our code differs from the standard ZEUS algorithm in two respects. First, we have 
implemented a version of the shearing-box boundary conditions. The model is then periodic 
in the y direction; the x boundaries are initial joined in a periodic fashion, but they are 
allowed to shear with respect to each other, becoming periodic again when t = nLy/{qQLx), 
n = 1, 2, . . .. A detailed description of the boundary conditions is given in Hawley, Gammie, 
& Balbus (1995). 

Second, we treat advection by the mean flow Vq = —qfixy separately from advection by 
the perturbed flow 6v = v — Vq. Mean-flow advection can be done by interpolation, using 
the algorithm described in Gammie (2001), which is similar to the FARGO scheme (Masset 
2000). This has the advantage that the timestep is not limited by the mean flow velocity (it 
is \6v\ rather than |t;| that enters the Courant condition). This permits the use of a timestep 
that is larger than the usual timestep by ~ L^/H if L^ ^ H. The shear-interpolation scheme 
also makes the algorithm more nearly translation-invariant in the x — y plane, thereby more 
nearly embodying an important symmetry of the underlying equations. 



2.1. Initial Conditions 

Without a specific model for the process that is injecting the vorticity, it is difficult 
to settle on a particular set of initial conditions, or to know how these initial conditions 
ought to vary when the size of the box is allowed to vary. Our choice of initial conditions is 
therefore somewhat arbitrary. We use a set of initial (incompressive) velocity perturbations 
drawn from a Gaussian random field. The amplitude of the perturbations is characterized 
by o" = (l^f/csp) . The power spectrum is |5wp ~ k~^^^, corresponding to the energy 
spectrum (E^ ~ k~^^^) of a two-dimensional Kolmogorov inverse turbulent cascade, with 
cutoffs at kmin = (V^) (27r/if ) and kmax = ^^'^kmin^- The surface density is not perturbed. 
These initial conditions correspond to a set of purely vortical perturbations. The parameters 
for our fiducial run are L^ = Ly = AH and a = 0.4. 



^We have compared our fiducial run to runs with a different range in fc, corresponding to vorticity injection 
either at scales ^ H or scales ~ O.IH. The results are qualitatively the same. 
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2.2. Code Verification 

Although our basic algorithm has already been tested (see Gammie 2001), we test the 
current version of our code by making a comparison with linear theory. Due to the underlying 
shear, small-amplitude perturbations in the shearing sheet are naturally decomposed in 
terms of shearing waves or shwaves (see Johnson & Gammie 2005 and references therein), 
Fourier components in the "co-shearing" frame. These have time- dependent wavenumber 
k{t) = kx{t)x + kyy, where kx{t) = k^o + q^kyt and k^o and ky are constant. The evolution of 
a single Fourier component can be calculated by integrating an ordinary differential equation 
for the amplitude of the shwave. For purely vortical (nonzero potential vorticity) or non- 
vortical perturbations, the evolution can be obtained analytically. The explicit expression 
for the amplitude of a vortical (incompressive) shwave is 

k^ 1 + r^ 

where k"^ = k^ + ky, r = qVtt + kxo/ky and a subscript on a quantity indicates its value 
at t = 0.^ The amplitude of a non-vortical (compressive) shwave satisfies the differential 
equation 

f^ + {cy + n')Svy, = o, (10) 

the solutions of which are parabolic cylinder functions. See Johnson & Gammie (2005) for 
further details on the shwave solutions. 

Figure 1 compares the numerical evolution of both vortical and compressive shwave am- 
plitudes with their analytic solutions. The initial shwave vector {kxo, ky) is {—16n/Lx, in/Ly) 
for the vortical shwave and {—Stt/Lx, 2t[/ Ly) for the compressive shwave. The other model 
parameters are the same as those in the fiducial run, except that L = 0.5H for the vortical- 
shwave evolution since kyH ^ 1 is required to prevent mixing between vortical and non- 
vortical shwaves near r = 0. The shwaves are well resolved until the radial wavelength 
Aa; = 4 X Ax, and the code is capable of tracking both potential-vorticity and compressive 
perturbations with high accuracy. 



3. Results 

The evolution of the potential vorticity in our fiducial run is shown in Figure 2. The 
snapshots are shown in lexicographic order beginning with the initial conditions, which have 



^This solution is valid at all times only for short-wavelength vortical perturbations (kH ^ 1). 
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equal positive and negative SC,- 

One of the most remarkable features of the fiducial run evolution is the appearance of 
comparatively-stable, long-lived vortices. These vortices have negative 6^ and are therefore 
dark in Figure 2. Similar vortices have been seen by Godon & Livio (1999, 2000), Li, Colgate, 
Wendroff, & Liska (2001) and Umurhan & Regev (2004). Cross sections of one of the vortices 
at the end of the run are shown in Figure 3. In our models the vortices are not associated 
with easily identifiable features in the surface density, since the perturbed vorticity is not 
large enough to require, through the equilibrium condition, an order unity increase in the 
local pressure. 

While the vortices are long-lived, they do decay. Figure 4 shows the evolution of the 
perturbed (noncircular) kinetic energy 

EK^ln6vl + 5vl) (11) 

in the fiducial run. Evidently the kinetic energy decays approximately as t~^/^ (which is 
remarkable in that, if the vortices would correspond to features in luminosity that decay 
as t~^'^, they could produce flicker noise; see Press 1978). Runs with half and twice the 
resolution decay in the same fashion, but if the resolution is reduced to 64^ the kinetic 
energy decays exponentially. Resolution of at least 128 zones per scale height appears to be 
required. 

What is even more remarkable is that the vortices are associated with an outward 
angular momentum flux, due to the driving of compressive motions by the vortices. Figure 5 
shows the evolution of the dimensionless angular momentum flux 

a = ^ / Tj6vjvydxdy (12) 

LxLyLlQCg J 

for models with a variety of resolutions. The data has been boxcar smoothed over an interval 
At = 10r2~^ to make the plot readable. Again, a resolution of at least 512^ appears to be 
required for a converged measurement of the shear stress. For the most highly resolved 
models a evolves like the kinetic energy, oc t~^/^. 

Compressibility is crucial for development of the angular momentum flux. We have 
demonstrated this in two ways. First, we have taken the fiducial run and decomposed 
the velocity field into a compressive and an incompressive part (i.e., into potential and 
solenoidal pieces in Fourier space) and measured the stress associated with each. For a 
set of snapshots taken from the last half of the fiducial run, the average total a = 0.0036; 
the incompressive component is a, = —0.0006; the compressive component is ac = 0.0032. 
The remaining alpha a^ = 0.00099 is in cross-correlations between the incompressive and 
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compressive pieces of the velocity field. As argued in Balbus (2000) and Balbus (2003), both 
incompressive trailing shwaves and incompressible turbulence transport angular momentum 
inward, whereas trailing compressive disturbances transport angular momentum outward. 
Our negative (positive) value for ai (ac) is consistent with this. 

Second, we have reduced the size of the model and reduced the amplitude of the initial 
perturbation so that it scales with the shear velocity at the edge of the model (constant 
"intensity" of the turbulence, in Umurhan and Regev's parlance). Thus the Mach number 
of the turbulence is reduced in proportion to the size of the box. We have compared four 
models, with L = (4, 2, 1, 0.5)if and a = (0.8, 0.4, 0.2, 0.1)cs. We would expect the lower 
Mach number models to have smaller-amplitude compressive velocity fields and therefore, 
consistent with the above results, smaller angular momentum fiux a. Averaging over the 
second half of the simulation, we find a = (0.0031, 0.0018, 7.2 x 10-^ -9.5 x 10"^). 

An additional confirmation of our overall picture can be seen in Figure 6, in which 
we show a snapshot of the velocity divergence superimposed on the potential vorticity for 
a medium-resolution (256^) version of the fiducial run.^ The position of the shocks with 
respect to the vortices in this figure is consistent with our interpretation that the former are 
generated by the latter. 

The smallest of our simulations {L = 0.5H) is nearly incompressible, but we continue 
to observe t'^^"^ decay (least squares fit power law is —0.49) at late times. The reason 
that we see decay while Umurhan & Regev (2004) do not may be that: (1) the remaining 
compressibility in our model causes added dissipation; (2) the numerical dissipation in our 
code is larger than that of Umurhan & Regev (2004); (3) the code used by Umurhan & Regev 
(2004) could somehow be aliasing power from trailing shwaves to leading shwaves (although 
they do explicitly discuss, and dismiss, this possibility). 

To highlight the dangers of aliasing for our finite-difference code, in Figure 7 we show 
the evolution of a vortical shwave amplitude at low resolution (64^), in units of r. We use the 
same parameters as those in our linear-theory test (Figure 1), for which the initial shwave 
vector corresponds to tq = —4. The initially-leading shwave swings into a trailing shwave, 
the radial wavelength is eventually lost near the grid scale, and due to ahasing the code 
picks up the evolution of the shwave again as a leading shwave. Repeating this test at higher 
resolutions indicates that successive swings from leading to trailing occur at an interval of 
r = N^/uy, where n^^ = 2 is the azimuthal wavenumber of the shwave. This is equivalent to 
kx{t) = 2tt/Ax. The decay of the successive linear solutions with time is due to numerical 



■^At higher resolutions, shocks are generated earher in the simulation from smaller vortices, and it is more 
difficult to see the effect we are describing due to the random nature of the vortices at this early stage. 



diffusion. 

Figure 7 suggests that it is easier to inject power into the simulation due to ahasing 
rather than to remove power due to numerical diffusion. We do not believe, however, that 
aliasing is affecting our high-resolution results. In addition, if we assume that the flow in our 
simulations can be modeled as two-dimensional Kolmogorov turbulence, then Svrms ~ A^^'^, 
where Svrms is the rms velocity variation across a scale A. The velocity due to the mean shear 
at these scales is Svghear ~ gf^A, and 5vrms/Sv shear ~ A"^/^. The velocities at the smallest 
scales are thus dominated by turbulence rather than by the mean shear. This conclusion is 
supported by the convergence of our numerical results at high resolution. 

Our model contains two additional numerical parameters: the size L and the initial 
turbulence amplitude a. Figure 8 shows the evolution of a for several values of a. Evidently 
for small enough values of a the a amplitude is reduced, but for near-sonic initial Mach 
numbers the a amplitude saturates (or at least the dependence on a is greatly weakened). 
Figure 9 shows the evolution for several values of L but the same initial a and the identical 
initial power spectrum. For large enough L the shear stress appears to be independent of L. 

Finally, we have studied the autocorrelation function of the potential vorticity as a means 
of characterizing structure inside the flow. Figure 10 shows the autocorrelation function 
measured in the flducial model and in an otherwise identical model with L = 8H. Evidently 
the potential vorticity is correlated over about one-half a scale height in radius, independent 
of the size of the model. This supports the idea that compressive effects limit the size of 
the vortices, since the shear flow becomes supersonic across a vortex of size ~ H (Barge & 
Sommeria 1995; Li, Colgate, Wendroff, & Liska 2001). 



4. Conclusion 

The presence of long-lived vortices in weakly-ionized disks may be an integral part of 
the angular momentum transport mechanism in these systems. The key result we have 
shown here is that compressibility of the flow is an extremely important factor in providing 
a signiflcant, positively-correlated average shear stress with its associated outward transport 
of angular momentum. Previous results using a local model have assumed incompressible 
flow and either report no angular momentum transport (Umurhan &; Regev 2004) or report 
a value {a ~ 10~^, Barranco & Marcus 2005) that is two orders of magnitude lower than 
what we find when we include the effects of compressibility. Global simulations (Godon & 
Livio 1999, 2000; Li, Colgate, Wendroff, & Liska 2001) have a difficult time accessing the high 
resolution that we have shown is required for a significant shear stress due to compressibility. 
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Our work leaves open the key question of what happens in three dimensions. Our 
vortices, which have radial and azimuthal extent < if, are inherently three-dimensional. 
Three-dimensional vortices are susceptible to the elliptical instability (Kerswell 2002) and 
are likely to be destroyed on a dynamical timescale. The fact that vortices persist in our two- 
dimensional simulations and not in the local (three-dimensional) shearing-box calculations 
of Balbus et al. (1996) is likely due to dimensionality. The recent numerical results of Bar- 
ranco & Marcus (2005) indicate that vortices near the disk midplane are quickly destroyed, 
whereas vortices survive if they are a couple of scale heights away from the midplane. Strong 
vertical stratification away from the midplane may enforce two-dimensional fiow and allow 
the vortices that we consider here to survive. 

The initial conditions in Barranco & Marcus (2005) are analytic solutions for two- 
dimensional vortices that are stacked into a three-dimensional column. The stable, off- 
midplane vortices apparently arise due to the breaking of internal gravity waves generated by 
the midplane vortices before they become unstable. There is also an unidentified instability 
that breaks a single off-midplane vortex into several vortices. These simulations leave open 
the question of whether stable off-midplane vortices can be generated from a random set of 
initial vorticity perturbations rather than the special vortex solutions that are imposed. 

Our work also leaves open the key question of what generates the initial vorticity. One 
possibility is that material builds up at particular radii in the disk, resulting in a global 
instability (e.g. Papaloizou & Pringle 1984, 1985) and a breakdown of the flow into vortices 
(Li, Colgate, Wendroff, & Liska 2001). A second possibility for vortex generation in variable 
systems is that the MHD turbulence, which likely operates during an outburst but decays 
as the disk cools (Gammie & Menou 1998), leaves behind some residual vorticity. The via- 
bility of such a mechanism could be tested with non-ideal MHD simulations such as those 
of Fleming & Stone (2003) and Sano & Stone (2003). A third possibility is that differential 
illumination of the disk somehow produces vorticity. Since the temperature of most cir- 
cumstellar disks is controlled by stellar illumination, small variations in illumination could 
produce hot and cold spots in the disk that interact to produce vortices. A fourth possibility 
is the generation of vorticity via baroclinic instability, which is likely to operate in disks 
whose vertical stratiflcation is close to adiabatic (Knobloch & Spruit 1986). The nonlinear 
outcome of this instability in planetary atmospheres is the formation of vortices, although it 
is far from clear that the same outcome will occur in disks. Finally, we note that a residual 
amount of vorticity can be generated from flnite-amplitude compressive perturbations. We 
have performed a series of runs with zero initial vorticity and perturbation wavelengths on 
the order of the scale height, and the results are qualitatively similar to Figure 5 with the 
shear stress reduced by nearly two orders of magnitude. 
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Fig. 1. — Evolution of velocity amplitudes for a vortical (left) and non- vortical (right) shwave. 
The heavy line is the analytic result, and the light lines are numerical results with (in order 
of increasing accuracy) N^j. = Ny = 32, 64, 128 and 256. 
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Fig. 2. — Panels show the evolution of the potential vorticity in the fiducial run. The 
size is 4H x 4H and the numerical resolution is 1024^. The initial conditions are shown 
in the upper left corner, and the other frames follow in lexicographic order at intervals 
of 22.2^2"^. Dark shades (blue and black in electronic edition) indicate potential vorticity 
smaller than Q/{2'Eq); light shades (yellow and red in electronic edition) indicate positive 
potential vorticity perturbations. Evidently only the "anticyclonic" (negative potential vor- 
ticity perturbation) vortices survive. Each vortex sheds sound waves, which steepen into 
trailing shocks. 
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Fig. 3. — Radial (left) and azimutlial (right) slice of a vortex at the end of the fiducial 
run. The heavy line shows the potential vorticity, the light line shows the magnitude of the 
velocity and the dotted line shows the surface density. 
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Fig. 4. — Evolution of kinetic energy in time for the fiducial run, on a log-log scale. The 
solid line shows a t~^/^ decay for comparison purposes. 
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Fig. 5. — Evolution of the shear stress a in the fiducial run and a set of runs at lower 
resolutions. 
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Fig. 6. — Snapshot of the velocity divergence superimposed on the potential vorticity in a 
medium- resolution (256^) version of the fiducial run. The thin (red in electronic edition) 
contours indicate negative divergence and are associated with shocks. The thick (blue in 
electronic edition) contours indicate negative potential vorticity. 
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Fig. 7. — Evolution of a vortical shwave amplitude in a low-resolution (64^) run, in units of 
r. The initial shwave vector {kxo,ky) is {—167f/Lx,4:7T/Ly), corresponding to tq = —4. The 
interval between successive peaks (a numerical effect due to aliasing) is r = N^/uy, where 
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Uy = 2 is the azimuthal wavenumber. 
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Fig. 8. — Evolution of the shear stress a in a set of runs at with varying initial a. Apparently 
for low values of cr the shear stress is reduced, but for initial Mach number near 1 the stress 
saturates. All runs have L = 4H. 
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Fig. 9. — Evolution of the shear stress a in a set of runs at with varying initial L, but the 
same initial Mach number a . 
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Fig. 10. — Autocorrelation function of the potential vorticity ^ for the fiducial model with 
L = AH (left) and a model with L = 8H (right). 



